package com.asa.numanaly;
/**
 * 逼近特征值
 * 做到这里，算法的难度越来越高，我的耐心也就
 * 所以实现最好用的，不会全实现了
 * @author Administrator
 *
 */
public class I {

	/**
	 * Householder算法
	 * 将对称三角矩阵转变成对角型的相似阵
	 * 注意这个必须是对三角矩阵，这个类里面还实现了一个针对任意矩阵的
	 * @param a
	 */
	public static double[][] tezhengzhi5(double[][] a) {
		double[][][] copya = new double[a.length][a.length][a.length];
		for (int i = 0; i < a.length; i++) {
			for (int j = 0; j < a.length; j++) {
				copya[0][i][j] = a[i][j];
			}
		}
		
		double[] v = new double[a.length];
		double[] u = new double[a.length];
		double[] z = new double[a.length];

		for (int k = 0; k < a.length-2; k++) {
			
			
			
			for (int i = 0; i < a.length; i++) {
				for (int j = 0; j < a.length; j++) {
					copya[k+1][i][j] = copya[k][i][j];
				}
			}
			
			
			
			
			
			
			
			double q = 0;
			for (int j = k+1; j < a.length; j++) {
				q += copya[k][j][k]*copya[k][j][k];
			}
			
			
			
			
			double alph;
			if (copya[k][k+1][k]==0) {
				alph = -Math.pow(q, 0.5);
			}else {
				alph = -Math.pow(q, 0.5)*copya[k][k+1][k]/Math.abs(copya[k][k+1][k]);
			}
			
			
			double RSQ = alph*(alph -copya[k][k+1][k]);
			
			
			v[k] = 0;
			v[k+1] = copya[k][k+1][k]-alph;
			
			for (int j = k+2; j < a.length; j++) {
				v[j] = copya[k][j][k];
			}
			
			
			
			
			for (int j = k; j < a.length; j++) {
				u[j] = 0;
				for (int i = k+1; i < a.length; i++) {
					u[j] = u[j]+copya[k][j][i]*v[i];
				}
			
				u[j] = u[j]/RSQ;
			}
			
			
			
			double PROD = 0;
			
			for (int i = k+1; i < u.length; i++) {
				PROD += v[i]*u[i];
			}
			
			
			
			for (int j = k; j < z.length; j++) {
				z[j] = u[j] - (PROD/(2*RSQ))*v[j];
				
			}
			
			
			for (int l = k; l < z.length; l++) {
				for (int j = l+1; j < z.length; j++) {
					
					copya[k+1][j][l] = copya[k][j][l]-v[l]*z[j] -v[j]*z[l];
					copya[k+1][l][j] = copya[k+1][j][l];
				}
				
				copya[k+1][l][l] = copya[k][l][l] -2*v[l]*z[l];
			}
			
			
			copya[k+1][a.length-1][a.length-1] = copya[k][a.length-1][a.length-1] - 2*v[a.length-1]*z[a.length-1]; 
			
			
			for (int j = k+2; j < z.length; j++) {
				copya[k+1][k][j] = 0;
				copya[k+1][j][k]=0;
			}
			
			copya[k+1][k+1][k] = copya[k][k+1][k] - v[k+1]*z[k];
			copya[k+1][k][k+1] = copya[k+1][k+1][k];

			
			
		}
		
		for (int i = 0; i < a.length; i++) {
			for (int j = 0; j < a.length; j++) {
				System.out.print(copya[copya.length-2][i][j]+"           ");
				
			}
			System.out.println();
		}
		return copya[copya.length-2];
		
		
		
	}
	
	
	
	
	/**
	 * Householder算法改进
	 * 这个能将非对称矩阵对三角化
	 * @param a
	 */
	public static double[][] tezhengzhi5_gai(double[][] a) {
		double[][][] copya = new double[a.length][a.length][a.length];
		for (int i = 0; i < a.length; i++) {
			for (int j = 0; j < a.length; j++) {
				copya[0][i][j] = a[i][j];
			}
		}
		
		double[] v = new double[a.length];
		double[] u = new double[a.length];
		double[] z = new double[a.length];
		double[] y = new double[a.length];

		for (int k = 0; k < a.length-2; k++) {
			
			
			
			for (int i = 0; i < a.length; i++) {
				for (int j = 0; j < a.length; j++) {
					copya[k+1][i][j] = copya[k][i][j];
				}
			}
			
			
			
			
			
			
			
			double q = 0;
			for (int j = k+1; j < a.length; j++) {
				q += copya[k][j][k]*copya[k][j][k];
			}
			
			
			
			
			double alph;
			if (copya[k][k+1][k]==0) {
				alph = -Math.pow(q, 0.5);
			}else {
				alph = -Math.pow(q, 0.5)*copya[k][k+1][k]/Math.abs(copya[k][k+1][k]);
			}
			
			
			double RSQ = alph*(alph -copya[k][k+1][k]);
			
			
			v[k] = 0;
			v[k+1] = copya[k][k+1][k]-alph;
			
			for (int j = k+2; j < a.length; j++) {
				v[j] = copya[k][j][k];
			}
			
			
			
			
			for (int j = 0; j < z.length; j++) {
				
				u[j] = 0;
				for (int i = k+1; i < a.length; i++) {
					u[j] = u[j]+copya[k][j][i]*v[i];
				}
			
				u[j] = u[j]/RSQ;
				
				
				y[j] = 0;
				for (int i = k+1; i < a.length; i++) {
					y[j] = y[j]+copya[k][i][j]*v[i];
				}
			
				y[j] = y[j]/RSQ;
				
				
			}
			double PROD = 0;
			for (int i = k+1; i < u.length; i++) {
				PROD += v[i]*u[i];
			}
			
			
			for (int j = 0; j < z.length; j++) {
				z[j] = u[j] - (PROD/RSQ)*v[j];
				
			}
			
			for (int l = k+1; l < z.length; l++) {
				for (int j = 0; j <= k; j++) {
					
					copya[k+1][j][l] = copya[k][j][l]-z[j]*v[l];
					copya[k+1][l][j] = copya[k][l][j]-y[j]*v[l];
				}
				for (int j = k+1; j < y.length; j++) {
					copya[k+1][j][l] = copya[k][j][l] -v[l]*z[j]-v[j]*y[l];

				}
			}
			
			
			

			
		}
		
		for (int i = 0; i < a.length; i++) {
			for (int j = 0; j < a.length; j++) {
				System.out.print(copya[copya.length-2][i][j]+"           ");
				
			}
			System.out.println();
		}
		
		return copya[copya.length-2];

		
	}
	
	
	/**
	 * QR算法	这个方法还可以直接得出特征向量，但要参考文档【Wil2】								未测试其正确性
	 * 20个步骤，头皮发麻
	 * 这个传入的必须是对称三角矩阵
	 * @param a
	 * @return
	 */
	public static double[][] tezhengzhi6(double[] a,double[] b,int M,double TOL) {
		int k = 0;
		double SHIFT=0;
		int n = a.length-1;
		
		double[][] copya = new double[M+1][a.length];
		double[][] copyb = new double[M+1][a.length];

		
		for (int i = 0; i < a.length; i++) {
			copya[0][i] = a[i];
		}
		for (int i = 0; i < b.length; i++) {
			copyb[0][i] = b[i];
		}
		
		
		
		
		
		while (k<=M) {
			
			if (Math.abs(copyb[k][n])<TOL) {
				double tzz = copya[k][n] +SHIFT;
				System.out.println(tzz);
				n=n-1;
			}
			
			
			if (Math.abs(copyb[k][2])<TOL) {
				double tzz = copya[k][1] +SHIFT;
				System.out.println(tzz);
				n = n-1;
				copya[k][1]=copya[k][2];
				for (int j = 2; j <= n; j++) {
					copya[k][j] = copya[k][j+1];
					copyb[k][j] = copyb[k][j+1];
				}
				
			}
			
			
			if (n==0) {
				return null;
			}
			
			if (n==1) {
				double tzz = copya[k][1]+SHIFT;
				System.out.println(tzz);
				return null;
			}
			
			
			for (int j = 3; j <= n-1; j++) {
				if (Math.abs(copyb[k][j])<TOL) {
					for (int i = 1; i <= j-1; i++) {
						System.out.print(copya[k][i]+"     ");
						
					}
					for (int i = 2; i <= j-1; i++) {
						System.out.print(copyb[k][i]+"     ");
						
					}
					System.out.println(" and ");
					for (int i = j; i <= n; i++) {
						System.out.print(copya[k][i]+"     ");
						
					}
					for (int i = j; i <= n; i++) {
						System.out.print(copyb[k][i]+"     ");
						
					}
					
					return null;

				}
			}
			
			
			//计算位移
			double bb = -(copya[k][n-1]+copya[k][n]);
			
			double cc = copya[k][n]*copya[k][n-1] -copyb[k][n]*copyb[k][n];
			double dd = Math.pow(bb*bb-4*cc,0.5);
			
			
			double u1;
			double u2;
			if (bb>0) {
				u1 = -2.0*cc/(bb+dd);
				u2 = -(bb+dd)/2;
		
			}else {
				u1 = (dd-bb)/2;
				u2 = 2.0*cc/(dd-bb);

			}
			
			
			
			if (n==2) {
				double tzz1 = u1+SHIFT;
				double tzz2 = u2+SHIFT;
				System.out.println(tzz1+"     "+tzz2);
				return null;
			}
			
			
			double s;
			
			if (Math.abs(u1-copya[k][n])>Math.abs(u2-copya[k][n])) {
				s = u2;
			}else {
				s = u1;
			}
			
			
			SHIFT = SHIFT+s;
			
			double[] c = new double[n+1];

			double[] d = new double[n+1];
			double[] r = new double[n+1];

			double[] sarr = new double[n+1];
			double[] q = new double[n+1];

			double[] x = new double[n+1];
			double[] y = new double[n+1];
			double[] z = new double[n+1];

			for (int j = 1; j <= n; j++) {
				d[j] = copya[k][j]-s;
			}
			
			x[1] = d[1];
			y[1] = b[2];

			
			for (int j = 2; j <= n; j++) {
				z[j-1]= Math.pow(x[j-1]*x[j-1]+copyb[k][j]*copyb[k][j], 0.5);
				
				
				c[j] = x[j-1]/z[j-1];
				sarr[j] = copyb[j-1][j]/z[j-1];
				q[j-1] = c[j]*y[j-1] + sarr[j]*d[j];
				x[j] = -sarr[j]*y[j-1]+c[j]*d[j];
				
				if (j!=n) {
					r[j-1] = sarr[j]*copyb[k][j+1];
					y[j] = c[j]*copyb[k][j+1];
				}
				
			}
			
			
			
			z[n] = x[n];
			copya[k+1][1] = sarr[2]*q[1]+c[2]*z[1];
			copyb[k+1][2] = sarr[2]*z[2];
			
			
			for (int j = 2; j < n-1; j++) {
				copya[k+1][j] = sarr[j+1]*q[j]+c[j]*c[j+1]*z[j];
				copyb[k+1][j+1] = sarr[j+1]*z[j+1];
			}
			copya[k+1][n] = c[n]*z[n];

			
			k++;
			
			
			
			
			
			
		}
		
		
		
		
		
		
		return null;
		
		
	}
	
	
	
	
	
	
	
	public static void main(String[] args) {
		
		double[][] a = 	{{4,1,-2,2},
						{1,2,0,1},
						{-2,0,3,-2},
						{2,1,-2,-1}};
		
		
		tezhengzhi5(a);
		System.out.println("=========================================");
		tezhengzhi5_gai(a);
	}
	
	
	
	
	
	
	
	
	
	
	
	
	
	
}
